function [P0,Y0,err,P] = new2dim(F,J,P0,delta,epsilon,max1)
% [P0,F0,err] = new2dim(F,J,P0,delta,epsilon,max1)
% [P0,F0,err,P] = new2dim(F,J,P0,delta,epsilon,max1)
% Newton's point iteration for higher dimensions.
% F is the vector function, input.
% J is the Jacobian matrix for F, input.
% P0 is the starting vector, input.
% delta is the tolerance for P0, input.
% epsilon is the tolerance for Y0, input.
% max1 is the maximum number of iterations, input.
% P0 is the vector solution, output.
% Y0 is the function vector, output.
% err is the error estimate for P0, output.
% P  is the matrix of iterations, output.
P = P0;
Y0 = feval(F,P0);
for k=1:max1,
  dF = feval(J,P0);
  if det(dF) == 0,
    dP = [0 0];
  else
    dP = (dF\Y0)';
  end
  P1 = P0 - dP;
  Y1 = feval(F,P1);
  err = norm(dP);
  relerr = err/(norm(P1)+eps);
  P0 = P1;
  Y0 = Y1;
  P = [P;P1];
  if (err<delta)|(relerr<delta)|(abs(Y1)<epsilon), break, end
end
Y0 = Y0';
